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LONG-TERM  GOALS 

The  long-term  goal  of  this  project  is  to  improve  the  existing  US  Navy  mesoscale  atmospheric  forecast 
model,  the  Coupled  Ocean/ Atmosphere  Mesoscale  Prediction  System  (COAMPS®1)  and  maintain  it  at 
the  leading  edge  of  the  Numerical  Weather  Prediction  (NWP)  models,  while  at  the  same  time  continue 
developing  the  next  generation  model.  In  the  near  future,  a  change  of  paradigm  to  unified  approach  in 
NWP  will  necessitate  a  novel  approach  in  modeling. 

OBJECTIVES 

There  are  two  main  objectives,  based  on  the  targets.  First,  to  continue  development  of  the  model  based 
on  the  Spectral  Element  (SE)  methods  and  evaluate  detailed  simulations  in  two-dimensions  (2D)  and 
preliminary  tests  in  three  dimensions  (3D).  In  addition,  improvements  to  the  numerical  temporal 
integration  scheme  will  be  sought.  Second,  to  improve  the  current  generation  version  of  CO  AMPS  by 
implementing  state-of-the-science  numerical  methods  such  as  Weighted  Essentially  Non-Oscillatory 
(WENO)  based  methods  for  tracer  transport  and  semi-lagrangian  time  stepping  for  improved  model 
efficiency.  The  overarching  objective  is  to  develop  a  dynamical  core  suitable  for  the  unified  approach 
in  support  of  NWP. 

APPROACH 

In  our  two-pronged  approach  we  will  leverage  our  efforts  to  advance  both  the  SE  model  and 
COAMPS. 

1.  SE  Model 

We  will  follow  a  systematical  plan  to  develop  and  test  the  new  model  based  on  the  SE  methods. 
Many  aspects  of  using  the  SE  model  in  the  atmospheric  science  are  new  so  consistency  of  the 
numerical  results  compared  to  known  solution  is  of  vital  importance.  The  SE  model  will  be 
evaluated  by  methodically  mapping  out  the  input  parameter  space  of  element  number  (h)  and 
polynomial  order  (p),  which  combined  determine  the  spatial  resolution.  We  will  use  two 
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idealized  two-dimensional  test  cases  to  assess  the  sensitivity  of  the  SE  model  to  the  setup 
parameters.  We  will  also  start  testing  the  code  in  three  dimensions.  In  addition,  a  few  new  time- 
integration  schemes  with  optimal  stability  properties  will  be  tested  as  a  possible  alternative  to 
the  semi-implicit  second  order  backward  difference  (BDF2)  method. 

2.  COAMPS 

Atmospheric  models  require  numerical  methods  that  can  accurately  represent  the  transport  of 
tracers  with  steep  gradients,  such  as  those  that  occur  at  cloud  boundaries  or  the  edges  of 
chemical  plumes.  In  atmospheric  sciences,  the  most  widely  used  numerical  techniques  for  this 
type  of  problem  are  flux-corrected  transport  or  closely  related  flux-limiter  methods.  The 
limiters  are  typically  designed  to  prevent  the  development  of  new  extrema  in  the  concentration 
field.  WENO  methods  are  widely  used  in  many  disciplines,  but  scarcely  been  tested  in  atmospheric 
applications.  Our  approach  has  been  to  carefully  test  the  WENO-based  scheme  for  a  variety  of 
idealized  and  real-data  scenarios. 

Key  personnel:  James  D.  Doyle  (Naval  Research  Laboratory),  P.  Alex  Reinecke  (NRL),  Sasa  Gabersek 
(NRL),  Dale  R.  Durran  (University  of  Washington),  Peter  N.  Blossey  (UW),  Francis  X.  Giraldo  (Naval 
Postgraduate  School),  Jim.  F.  Kelly  (NPS). 

WORK  COMPLETED 

We  have  completed  an  in-depth  analysis  of  the  sensitivity  of  the  simulation  results  to  the  setup 
parameters  for  two  idealized  two-dimensional  test  cases.  The  h-p  parameter  space  was  systematically 
mapped  out  resulting  in  91  pairs  of  setup  parameters  h  and  p.  The  average  spatial  resolution,  defined  as 
the  average  nodal  spacing  within  the  element,  varied  from  200  m  to  10  km  in  horizontal  and  40  m  to 
1500  m  in  vertical  direction.  The  first  test  represents  a  flow  over  a  small  hill,  resulting  in  a  linear, 
hydrostatic  mountain  wave  for  which  an  analytic  solution  is  known  and  simulation  results  are 
quantitatively  assessed.  The  second  test  is  a  simulation  of  an  idealized  mid-latitude  squall  line, 
triggered  by  a  thermal  bubble  at  the  start  of  the  simulation.  The  storm  evolution  and  characteristic 
quantities,  such  as  precipitation  accumulation,  are  evaluated.  We  have  also  completed  a  preliminary 
series  of  three-dimensional  simulations  of  an  idealized  flow  over  an  isolated  mountain.  Four 
alternative  implicit-explicit  time  integration  schemes  with  desirable  overall  stability  properties  were 
tested  on  a  few  idealized  two-dimensional  cases  and  compared  to  the  existing  scheme  (BDF2)  in  terms 
of  maximum  stable  time  step  and  elapsed  wallclock  time. 

The  WENO  based  selective  monotonic  advection  (SMA)  scheme  (Blossey  and  Durran,  2007)  has  been 
fully  implemented  in  COAMPS.  In  addition,  a  semi-Lagrangian  option  (Skamarock,  2006;  Blossey 
and  Durran,  2007)  has  also  been  implemented  in  order  to  offset  the  extra  computational  costs  required 
for  the  SMA  scheme.  Both  the  semi-Lagrangian  and  Eulerian  versions  of  the  SMA  scheme  were 
tested  and  compared  with  the  currently  available  2nd-  and  4th-order  finite  difference  advection 
schemes  for  a  variety  of  idealized  and  real  data  test  problems. 

RESULTS 

The  results  of  the  SE  model  for  a  linear,  hydrostatic  mountain  wave  suggest  that  the  model  is 
computationally  more  expensive  at  the  same  resolution  compared  to  a  finite-difference  model.  If  an 
apriori  acceptable  value  of  error  is  known,  the  SE  model  is  less  expensive  because  it  achieves  the 
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desired  accuracy  at  a  coarser  resolution.  Refining  the  spatial  resolution  by  increasing  the  polynomial 
order,  number  of  elements,  or  both,  becomes  computationally  expensive  and  is  no  longer  justified  as 
soon  as  the  nominal  grid  spacing  is  adequate  to  resolve  the  simulated  phenomenon. 

The  squall  line  simulations  confirm  the  robustness  of  the  SE  model  across  the  parameter  space,  when 
the  storm  is  adequately  resolved.  With  no  analytic  solution  available,  we  compared  the  storm  evolution 
and  some  selected  integrated  quantities.  When  the  average  nodal  spacing  is  adequate  to  describe  the 
physical  phenomenon,  the  results  over  the  parameter  space  look  very  similar  (Fig  1),  with  similar 
conclusions  for  other  values  of  the  nominal  spacing. 


Figure  1:  Vertical  cross  sections  of  the  squall  line  evolution  at  t=6000  s  with  the  nominal  spacing 
Ax=l  km.  The  polynomial  order  varies  as  follows:  p=4  (top  left  panel),  p=6  (top  right  panel),  p=8 
(bottom  left  panel),  p=10  (bottom  right  panel).  Filled  contours  represent  equivalent  potential 
temperature  perturbation  (contour  interval  3  K),  with  red  (blue)  shading  representing  positive 
(negative)  values.  The  interval  centered  around  zero  (1-3,  3])  is  omitted.  The  cloud  mixing  ratio 
(qc=10's)  in  thick  black  line  represent  the  cloud  outline.  The  bottom  portion  of  each  panel  shows 

rain  water  accumulation  as  a  function  of  distance. 


Unlike  for  the  dry  simulations,  refining  the  nominal  resolution  does  not  result  in  convergence  of 
solutions.  A  noticeable  trend  in  some  integrated  quantities,  such  as  a  decrease  of  rainfall  accumulation 
(Fig  2)  and  maximum  rain  rate  with  the  increased  resolution  and  an  increase  of  maximum  vertical 
velocity  can  be  attributed  to  nonlinearity  of  moist  processes.  As  the  grid  spacing  gets  finer,  the  scale  of 
instabilities  is  refined  as  well.  Instead  of  a  large  updraft,  there  are  many  smaller  cells  with  more  intense 
vertical  velocities.  Their  overall  combined  effect,  in  terms  of  the  mass  flux  or  transport  of  moisture 
from  the  surface  layer  into  the  convective  cell,  and  subsequent  cloud  processes  (e.g.  water  vapor 
condensation  with  the  latent  heat  release),  is  smaller  compared  to  a  larger  updraft,  resulting  in  reduced 
precipitation  rates  and  consequently  less  accumulated  precipitation. 
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Figure  2:  Total  rainfall  accumulation  (mm)  in  6  hours  averaged  over  the  entire  domain  as  a 
function  of  the  polynomial  order  (p)  and  number  of  elements  in  horizontal  (h)  for  the  squall  line 
simulations.  Curved  color  lines  represent  constant  nominal  horizontal  resolution  (Ax,  or  constant 
number  of  nodal  points  in  the  horizontal  direction):  green  (Ax=3  km),  yellow  (Ax=2  km),  orange 

(Ax=l  km)  and  red  (Ax=0.5  km). 

To  provide  further  evidence,  we  performed  a  horizontal  spectral  analysis  of  the  vertical  velocity, 
averaged  over  time  and  height  (Fig  3).  In  the  original  set  of  simulations  (Fig  3,  left  panel),  we  noticed 
that  with  the  increasing  resolution  the  peak  of  the  power  spectrum  moved  towards  shorter  scales.  We 
performed  an  additional  set  of  “dry”  simulations  with  no  moisture  and  the  resulting  power  spectra 
coincide  regardless  of  the  spatial  resolution  (Fig  3,  right  panel). 


Figure  3:  Power  spectra  for  simulations  with  p=8  and  varying  nominal  resolutions  Ax=3  km 
(thickest,  lightest  blue),  Ax=1.5  km  (thick,  light  blue),  Ax=0. 75  km  (thin,  dark  blue)  and  Ax=0.375 
km  (thinnest,  darkest  blue).  The  left  panel  is  for  the  control  squall  line  simulations  and  the  right 
panel  is  for  the  “dry”  squall  line  with  no  moisture.  The  spectra  are  averaged  over  height  (0-12  km, 
with  0.5  km  increment)  and  time  (0-4  h,  with  300  s  increment). 
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The  SE  model  has  been  expanded  and  is  currently  being  tested  in  three  dimensions  (e.g.,  flow  over  an 
isolated  mountain  is  shown  in  Fig.  4).  In  addition,  preliminary  tests  of  scaling  on  massively  parallel 
computational  platforms  are  also  being  evaluated. 


Figure  4:  Linear,  inviscid  hydrostatic  flow  over  an  idealized  3D  mountain  (circular  terrain  contours 
are  visible  in  the  x-y  plane).  Isosurfaces  of  vertical  velocity  perturbation  with  red  (blue)  color 
representing  positive  (negative)  value  of  0.001  m  s-1.  The  full  horizontal  extent  of  the  domain  with 
element  structure  in  the  background  (left  panel)  and  a  detailed  cross-section  at  the  center  of  the 
domain  showing  the  nodal  points  within  a  few  elements. 

Four  additional  time  integration  schemes  were  tested  using  an  idealized  2D  test  case  (Schaer 
mountain).  Some  candidate  schemes  allowed  for  approximately  double  the  time  step  permitted  by  the 
currently  used  BDF2  scheme,  but  the  overall  time  to  finish  the  simulation  was  in  a  couple  of  cases 
much  longer  (bdf23n  and  am2sax32).  It  is  interesting  to  note  that  the  overall  fastest  execution  time  was 
achieved  by  the  mcnax32  scheme  when  using  the  same  time  step  as  bdf2.  Longer  time  steps  cause 
more  iterations  while  solving  a  system  of  linear  equations. 

Table  1:  Time  integration  schemes  with  the  maximum  stability  permissible  time  step  and  total 

execution  time  for  a  simulation  lasting  50000s. 


Scheme  name 

Maximum  stable  time  step  dt 

Walltime  to  execute 

bdf2  (default) 

dt  =  4  s 

2379  s 

bdf23n 

dt  =  7.6  s 

3142  s 

am2sax32 

dt  =  8.4  s 

3981  s 

mcnax32 

dt  =  7.6  s 

2120  s 

mcnax32 

dt  =  4  s 

1581  s 
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Figure  5  shows  results  from  a  series  of  two-dimensional  mid-latitude  squall  line  simulations  for  the 
WENO  based  selective  monotonic  advection  scheme  in  COAMPS.  Plotted  is  a  one-hour  average 
vertical  cross  section  of  simulated  cloud-liquid  water,  rain  water,  and  negative  potential  temperature 
perturbations.  While  differences  exist  between  the  4th-order  advection  scheme  (Fig.  5a)  and  the  SMA 
scheme  (Fig.  5b),  the  overall  qualitative  behavior  is  similar.  Both  the  4th-order  advection  and  the  SMA 
scheme  contain  a  compact  convective  updraft,  spreading  low-level  cold  pool  near  the  surface,  and  a 
broad  anvil  in  the  upper  troposphere.  Figures  5cd  show  results  from  a  similar  two-dimensional  squall¬ 
line  simulation  but  with  the  semi-lagrangian  scheme.  The  simulation  with  a  time  step  two-times  as 
long  as  the  control  (Fig.  5c)  is  qualitatively  similar  to  the  Eulerian  versions  of  the  code.  On  the  other 
hand,  as  the  length  of  the  time  step  increases,  intensity  of  the  squall  decreases  and  the  spreading  of  the 
anvil  is  limited  in  the  upper  troposphere  (Fig.  5d). 

One  drawback  of  the  SMA  scheme  is  the  extra  computational  overhead  required  to  maintain 
monotonicity.  For  example,  the  SMA  squall-line  simulation  requires  approximately  25%  more  cpu 
time  than  the  similar  4th-order  advection  solution.  Fortunately,  switch  from  the  Eulerian  framework 
to  the  semi-lagrangian  framework  can  help  minimize  the  extra  computational  burden.  The  run  time  for 
the  4dt  semi-lagrangian  simulation  is  only  12%  longer  than  the  4th-order  advection  solution.  It  is 
expected  further  improvements  in  timing  will  be  realized  for  the  semi-lagrangian  scheme  as  more 
scalar  variables  are  included  and  the  domain  expanded  to  three  dimensions. 
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Figure  5.  One-hour  average  vertical  cross  section  of  a  two-dimensional  squall  simulation  with  (a) 
fourth-order  advection,  (b)  selective  monotonic  advection,  (c)  selective  monotonic  advection  with  a 
semi-lagrangian  time  step  that  is  2-  times  the  time  step  of  (b),  and  (d)  same  as  (c)  except  the  time 
step  is  4-times  the  time  step  of  (b).  Plotted  is  the  mixing  ratio  for  cloud  liquid  water  (filled,  ci=0.1 
g/kg)  and  0.25  g/kg  rain  water  contour  (white)  as  well  as  negative  perturbation  potential 

temperature  (blue;  ci=2  K). 
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IMPACT/APPLICATIONS 


COAMPS  is  the  Navy’s  operational  mesoscale  NWP  system  and  is  recognized  as  the  key  model 
component  driving  a  variety  of  DoD  tactical  decision  aids.  Accurate  mesoscale  prediction  is 
considered  an  indispensable  capability  for  defense  and  civilian  applications.  Skillful  CO  AMPS 
predictions  at  resolutions  less  than  1  km  will  establish  new  capabilities  for  the  support  of  the 
warfighter  and  Sea  Power  21.  The  SE  model  could  become  the  unified  dynamical  core  for  both  global 
and  mesoscale  weather  forecasts  of  the  US  Navy.  The  design  of  the  model  and  its  structure  are  such 
that  it  will  take  a  full  advantage  of  available  computational  cores  as  they  become  available  in  large 
numbers  (-100000)  in  the  near  future. 

TRANSITIONS 

The  next  generation  COAMPS  system  will  transition  to  6.4  projects  within  PE  0603207N  (SPAWAR, 
PMW-120)  that  focus  on  the  transition  CO  AMPS  to  FNMOC.  The  improvements  to  the  COAMPS 
dynamical  core  have  been  transitioned  to  the  SPAWAR  6.4  project  and  subsequently  to  operations  as  a 
result  of  the  marked  improvement  in  the  geopotential  height  bias  statistics. 

RELATED  PROJECTS 

COAMPS  will  be  used  in  related  6.1  projects  within  PE  0601 153N  that  include  studies  of  air-ocean 
coupling,  boundary  layer  studies,  and  topographic  flows  and  in  related  6.2  projects  within  PE 
0602435N  that  focus  on  the  development  of  the  atmospheric  components  (QC,  analysis,  initialization, 
and  forecast  model)  of  COAMPS. 
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